DSMC  Simulations  in  Support  of  the  STS- 107  Accident 

Investigation 

Michael  A.  Gallis+ 

Sandia  National  Laboratories,  Albuquerque,  New  Mexico  87185-0826,  USA 

% 

Katie  A.  Boyles  &  Gerald  J.  LeBeau 

NASA  Johnson  Space  Center,  Houston,  Texas  77058-EG3,  USA 

Abstract.  Three-dimensional  Direct  Simulation  Monte  Carlo  simulations  of  Columbia  Shuttle  Orbiter  flight 
STS- 107  are  presented.  The  aim  of  this  work  is  to  determine  the  aerodynamic  and  heating  behavior  of  the 
Orbiter  during  aerobraking  maneuvers  and  to  provide  piecewise  integration  of  key  scenario  events  to  assess 
the  plausibility  of  the  candidate  failure  scenarios.  The  flight  of  the  Orbiter  is  examined  at  two  altitudes:  350 
kft  and  300  kft.  The  flow  field  around  the  Orbiter  and  the  heat  transfer  to  it  are  calculated  for  the  undamaged 
configuration.  The  flow  inside  the  wing  for  an  assumed  damage  to  the  leading  edge  in  the  form  of  a  10-inch 
hole  is  studied. 


INTRODUCTION 

On  January  16,  2003,  the  Space  Shuttle  Columbia  was  launched  from  Kennedy  Space  Center  to  begin  STS- 
107,  the  113th  mission  in  the  Space  Shuttle  Program  and  Columbia’s  28th  flight.  As  it  passed  just  northwest  of 
Hawaii  sixteen  days  later  on  Saturday,  February  1,  2003,  Columbia  initiated  a  maneuver  to  return  to  earth  for  a 
landing  at  Kennedy  Space  Center  in  Florida.  By  the  time  it  reached  the  US  mainland  and  continued  its  descent 
over  the  western  United  States,  people  on  the  ground  observed  bright  objects  coming  off  of  Columbia. 
Columbia  broke  apart  at  an  altitude  of  approximately  200,000  feet  near  Dallas,  Texas  at  approximately  8  AM 
CST,  just  15  minutes  before  its  scheduled  landing  at  Kennedy  Space  Center  in  Florida. 

The  conclusion  of  the  investigation  [1]  was  that  the  physical  cause  of  the  loss  of  Columbia  and  its  crew  was 
a  breach  in  the  Thermal  Protection  System.  At  81.7  seconds  after  the  launch,  a  piece  of  insulating  foam 
separated  from  the  left  bipod  ramp  section  of  the  External  Tank  and  struck  the  Orbiter,  causing  a  breach  on  the 
leading  edge  of  the  left  wing.  During  re-entry,  this  breach  in  the  Thermal  Protection  System  allowed 
superheated  air  to  penetrate  through  the  leading  edge  insulation  and  progressively  melt  the  aluminum  structure 
of  the  left  wing,  resulting  in  a  weakening  of  the  structure  until  increasing  aerodynamic  forces  caused  loss  of 
control,  failure  of  the  wing,  and  break-up  of  the  Orbiter.  The  investigation  of  this  scenario  involved  numerous 
CFD  [1]  and  Direct  Simulation  Monte  Carlo  (DSMC)  [1,2]  simulations  to  provide  the  heating  and  forces  on  the 
Orbiter  and,  in  particular,  on  the  leading  edge  of  the  left  wing.  In  this  paper,  representative  DSMC  calculations 
conducted  in  support  of  the  investigation  of  this  particular  failure  scenario  will  be  presented.  It  should  be  noted 
that  at  the  time  when  this  work  was  performed  the  conclusion  of  the  investigation  was  not  available.  The 
purpose  of  the  simulations  presented  herein  was  to  assist  the  work  of  the  investigation  committee,  especially  in 
an  area  where  forensic  and  physical  data  were  not  available. 

In  the  first  part  of  the  paper,  the  undamaged  full  3-D  geometry  is  studied  and  the  flow  field  characteristics 
analyzed.  In  the  second  part  of  the  paper,  boundary  conditions  are  extracted  from  the  full  body  simulations  to  be 
used  in  examining  the  wing  leading  edge  with  an  assumed  10-inch  diameter  hole.  It  is  interesting  to  note  that 
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after  this  analysis  was  completed,  an  experimental  reproduction  of  the  impact  on  the  wing  leading  edge 
indicated  that  the  size  of  the  hole  might  have  been  even  greater  than  10  inches  and  possibly  as  large  as  16 
inches. 


METHODOLOGY  USED 

For  the  purposes  of  this  work,  the  "DAC"  DSMC  [2]  implementation  of  LeBeau  [3]  was  used.  Since 
the  DSMC  algorithm  is  well  established  on  physical  arguments,  DSMC  implementations  differ  mostly  in  the 
flow  field  discretization  methods  and  the  way  they  represent  the  modeled  geometry.  The  DAC  software 
employs  a  two-level  embedded  Cartesian  grid  system.  Thus  the  computational  domain  is  a  rectangular  box, 
aligned  with  the  Cartesian  axes.  Embedded  within  the  flow  field  grid  is  the  surface  geometry.  DAC  represents 
the  surface  geometry  as  a  collection  of  unstructured  triangular  elements,  which  also  act  as  sampling  zones  for 
surface  properties. 

The  rectangular  bounding  box  for  the  computational  domain  is  specified  by  the  user,  as  is  the  discretization 
in  each  of  the  three  Cartesian  directions.  The  cells  created  by  this  uniform  Cartesian  grid  are  referred  to  as 
Level-I  Cartesian  cells.  Each  of  these  Level-I  Cartesian  cells  can  be  further  refined  by  its  own  embedded 
Cartesian  grid.  These  embedded  Cartesian  grids  form  Level-II  Cartesian  cells.  The  two-level  embedded 
Cartesian  grid  system  permits  variable  refinement  throughout  the  computational  domain,  which  is  essential  for 
meeting  the  local  mean-free-path  cell  size  requirement.  To  model  chemical  reactions,  Bird’s  Total  Collision 
Energy  model  (TCE)  model  was  used.  The  chemical  reaction  data  used  in  the  simulations  are  those  proposed  by 
Park  [4]. 


SHUTTLE  ORBITER  MODEL  DEVELOPMENT 

The  model  of  the  Shuttle  Orbiter  was  based  on  a  detailed  CAD  design  provided  by  Boeing  Huntington 
Beach.  From  this,  a  surface  grid  of  approximately  32,000  triangular  elements  was  constructed.  A  Cartesian  gas- 
phase  grid  of  more  than  5,000,000  elements  (after  the  adaptation)  was  used  for  the  computational  domain. 
Every  possible  effort  was  made  to  ensure  that  all  these  parameters  complied  with  Bird’s  criteria  [2]  for  a 
successful  DSMC  simulation.  The  time  step  was  such  that  molecules,  on  average,  did  not  move  more  than  a  cell 
per  time  step,  and  cells  were  smaller  than  the  local  mean  free  path.  For  the  gas-phase  grid  adaptation,  the 
“Level-II  cell”  procedure  (outlined  in  LeBeau  [3])  was  used.  The  adaptation  procedure  allowed  up  to  20  Level- 
II  cells  in  each  Level-I  cell. 

The  solution  was  run  fully  diffuse,  and  the  surface  wall  temperatures  used  for  the  simulations  were 
estimated  to  be  525  K  for  the  350-kft  case  and  950  K  for  the  300-kft  case.  A  chemistry  model  for  high- 
temperature  reacting  air  was  used  that  contained  six  molecular  species:  O^,  N2,  O,  N,  NO,  and  Ar.  The  mole 
fractions  of  these  species  were  obtained  from  the  1976  Standard  Atmosphere  Model  [5].  Number  density, 
velocity,  and  free  stream  temperature  were  obtained  from  STS- 107  flight  data. 

FLOW  FIELD  CALCULATIONS  OF  THE  UNDAMAGED  GEOMETRY 

For  the  purposes  of  the  investigation,  two  cases  were  examined.  The  first  one  was  in  the  rarefied  regime 
and  the  second  the  transitional  regime.  The  conditions  of  the  two  instances  are  given  in  Table  1.  For  these 
conditions  no  CFD  calculations  were  performed.  The  times  corresponding  to  the  two  cases  are  measured  in 
seconds  from  “Entry  Interface”  (El),  an  arbitrarily  determined  altitude  of  400  kft,  where  the  Orbiter  begins  to 
experience  the  effects  of  Earth’s  atmosphere.  Entry  Interface  for  STS- 107  occurred  at  8:44:09  AM  (EST)  on 
February  1st,  2003.  The  first  abnormal  indication  occurred  270  seconds  after  Entry  Interface,  i.e.  approximately 
70  seconds  after  the  second  instance  (300  kft)  occurred.  This  fact  was  known  to  neither  the  crew  nor  Mission 
Control  since  sensor  data  were  recorded  and  stored  onboard  the  Orbiter. 

Figures  1  and  2  present  the  average  temperature  (average  flow  energy)  profiles  around  the  Orbiter  for 
the  350-kft  and  the  300-kft  cases,  respectively.  The  flow  travels  from  right  to  left.  In  the  350-kft  case,  as  the 
first  free  stream  molecules  hit  the  surface  of  the  vehicle,  they  are  slowed  down,  and,  as  a  result,  a  large 
compression  region  is  formed.  This  compression  region  is  an  extended  Knudsen/boundary  layer.  Diffusion  of 
air  molecules  upstream  (due  to  the  density  gradient)  gives  rise  to  head-on  collisions  between  these  molecules 
and  the  free-stream-air  incoming  molecules.  These  highly  energetic  collisions  lead  to  an  increase  of  the 
temperature  of  the  flow  that  is  clearly  shown  as  the  extended  high-temperature  area.  The  thickness  of  this  area 
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is  of  the  order  of  ten  local  mean  free  paths.  As  the  flow  moves  closer  to  the  surface  of  the  vehicle,  the 
temperature  of  the  flow  decreases  since  for  the  simulation,  the  temperature  of  the  vehicle  was  kept  constant  for 
each  instant. 


TABLE  1.  Flight  conditions  for  the  DSMC  cases  examined  during  the  STS-107  investigation 
Properties_ 350  kft_ 300  kft 


Time  from  Entry  Interface  (s) 

91 

197 

Mach  Number 

25.1 

27.0 

Altitude  (km) 

106.7 

91.4 

Altitude  (ft) 

350,274 

300,003 

Angle  of  Attack  (deg) 

41 

40 

Knudsen  number 

0.02 

0.001 

In  Figure  2,  the  temperature  profile  is  shown  for  the  300-kft  case.  A  distinct  shock  layer  is  developed 
starting  from  the  nose  of  the  vehicle  and  expanding  along  it.  The  temperature  reaches  almost  15,000  K.  It  is 
worth  noting  that,  although  the  shock  layer  is  much  thinner  in  this  case  due  to  the  shorter  mean  free  path,  the 
temperatures  achieved  are  similar  to  those  in  the  previous  case.  The  higher  density  of  this  case  (by  almost  one 
order  of  magnitude)  resulted  in  a  heat  flux  greater  than  the  350-kft  case  by  almost  one  order  of  magnitude. 

In  the  wake  of  the  vehicle,  both  the  density  (not  shown  here)  and  the  temperature  drop  by  one  order  of 
magnitude  due  to  the  rapid  expansion.  More  details  about  these  flow  fields  can  be  found  in  Reference  [6]. 
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FIGURE  1.  Temperature  profile  (350-kft  case) 


FIGURE  2.  Temperature  profile  (300-kft  case) 


FLOW  THROUGH  A  10-INCH  HOLE  IN  THE  LEADING  EDGE 

In  the  second  part  of  the  analysis,  the  internal  flow  behind  the  Re-enforced  Carbon-Carbon  (RCC)  panels  7, 
8,  and  9  due  to  a  10-inch  hole  in  RCC  panel  8  was  simulated.  For  the  simulation  of  the  internal-flow  cases,  the 
inflow  boundary  conditions  were  extracted  from  the  previously  run  DSMC  external  flow  solutions.  Vents  were 
also  applied  to  each  side  of  the  geometry  with  a  total  area  equivalent  to  66  in2  to  emulate  the  venting  holes 
along  the  upper  surface  of  the  wing. 

A  heritage  CAD  definition  for  the  RCC  cavity  was  obtained  from  Jim  Greathouse  [7].  This  definition 
includes  wing  geometry  from  the  leading  edge  back  to  the  wing  spar  and  includes  a  definition  for  the  earmuffs 
between  panels.  The  geometry  used  in  this  case  was  developed  using  the  GridGen  grid  generation  tool.  After 
importing  the  CAD  into  GridGen,  a  10-inch  hole  was  generated  in  RCC  panel  8  at  a  location  of  X  =  1065 
inches,  Y  =  -219  inches,  and  Z  =  286  inches  in  the  Orbiter  coordinate  system.  Vent  holes  were  included  on  each 
side  of  the  geometry.  It  was  assumed  that  the  RCC  thickness  was  uniform  at  14  inch,  so  the  10-inch  hole  in  RCC 
panel  8  was  given  a  thickness  of  approximately  %  inch. 
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Since  DSMC  solutions  had  already  been  obtained  for  the  undamaged  cases,  it  was  felt  that  it  was  not 
necessary  for  a  full-length  shuttle  case  to  be  run  for  this  failure  scenario.  Instead,  an  external-flow  field  box-like 
’’geometry”  was  created  that  connected  to  the  hole  from  the  external  flow  side.  Great  care  was  taken  to  create  an 
external  flow  field  geometry  large  enough  so  that  the  geometry  of  the  hole  would  not  significantly  affect  the 
flow  at  the  upstream  boundary  but  small  enough  so  that  the  case  would  run  in  a  timely  manner.  The  external 
flow  field  geometry  was  created  to  connect  to  a  corresponding  portion  of  the  outer  RCC  surface  which,  in  turn, 
connected  to  the  surface  representing  the  thickness  of  the  hole. 

In  DAC,  the  surface  geometry  is  represented  as  a  collection  of  unstructured  triangular  elements  and  the 
entire  surface  geometry  must  be  created  to  be  ’’watertight”  and  possess  continuous  connectivity  among  triangles. 
In  this  case,  the  surface  geometry  contained  approximately  53,000  surface  triangles.  Each  individual  surface 
triangle  has  specific  boundary  conditions  associated  with  it.  The  triangles  representing  the  surfaces  of  the  wing 
were  given  the  ’’solid  wall”  boundary  condition,  which  means  that  no  molecules  are  produced  at  the  wall  and 
none  pass  through.  The  triangles  representing  the  vents  were  given  the  ’’outflow”  boundary  condition,  which 
means  that  any  molecules  that  hit  the  surface  will  disappear.  Finally,  the  triangles  representing  the  external  flow 
field  surface  were  given  the  ’’inflow”  boundary  condition.  An  inflow  boundary  creates  molecules  on  the  flow 
field  side  of  the  triangle,  but  any  molecules  that  hit  the  surface  will  disappear. 

The  computational  domain  was  specified  to  encompass  the  surface  geometry  created  in  this  case  and  the 
boundary  conditions  of  the  computational  domain  were  specified  as  vacuum.  For  the  300-kft  case,  the  Level-I 
discretization  employed  cells  of  0.0015  meters.  The  total  number  of  cells  after  the  adaptation  was  in  excess  of 
10  million  (grids  up  to  40  million  cells  were  also  used).  For  the  350-kft  case,  the  Level-I  cells  were  sized  at 
0.033  meters.  The  total  number  of  cells  after  the  adaptation  was  approximately  1  million. 

The  wall  temperature  was  assumed  to  be  300  K.  The  chemistry  model  used  was  the  same  as  for  the 
external-  flow  cases. 


Internal  Flow  Results 


Figure  3  displays  a  Y-slice  of  the  flow  field  number  density  for  the  internal  flow  case  at  350  kft.  In  this 
figure  the  streamlines  shown  depict  the  movement  of  the  flow  as  it  enters  the  hole  in  panel  8  and  begins 
recirculating  inside  the  cavity.  The  exact  location  of  these  panels  along  the  leading  edge  of  the  vehicle  is  given 
in  Reference  1.  In  this  simulation,  the  flow  is  allowed  to  escape  through  the  exits  at  panels  7  and  9;  however, 
the  total  area  provided  was  not  enough  for  the  flow  to  exit  and  as  a  result  the  recirculation  zone  shown  in  Figure 
3  is  formed.  A  result  of  the  formation  of  a  recirculation  zone  in  the  leading  edge  cavity  is  that  the  heat  transfer 
to  the  wing  spar  is  increased.  Another  undesirable  result  of  the  existence  of  flow  in  the  wing  cavity  is  to  weaken 
the  integrity  of  the  wing  spar  by  removing  layers  of  material  from  it. 
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FIGURE  3.  Number  density  and  streamlines  for  the 
350-kft  internal- flow  case 


FIGURE  4.  Reference  heating  for  the  350-kft 
internal-flow  case 


Figure  4  shows  the  normalized  (by  the  stagnation  point  heat  flux)  heating  distribution  in  the  RCC  cavity  for 
the  350-kft  case.  The  heating  values  shown  in  the  legends  on  the  right  side  of  the  figure  have  been  normalized 
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by  reference  heating  values  based  on  the  free-stream  properties.  From  these  figures,  it  is  evident  that  the  area  of 
highest  heating  due  to  the  breach  is  near  the  area  between  panels  8  and  9.  As  the  flow  hits  the  comer  of  the 
earmuff,  a  shock  is  formed  and  the  flow  is  seen  to  splash  onto  the  spar.  Elevated  heating  levels  are  also  seen 
between  panels  7  and  8. 
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FIGURE  5.  Number  density  and  streamlines  for  the 
300-kft  internal-flow  case 


FIGURE  6.  Reference  heating  for  the  300-kft 
internal-flow  case 


Figures  5  and  6  show  the  number  density  and  normalized  heating  profiles  for  the  300-kft  case.  In  harmony 
with  the  350-kft  case,  the  initial  free-stream  conditions  were  extracted  from  the  300-kft  simulation.  Figure  5 
shows  the  density  and  flow  streamlines.  In  a  similar  fashion  as  the  350-kft  case,  a  recirculation  region  is  formed 
inside  the  wing  cavity.  The  normalized  surface  heating  is  given  in  Figure  6.  The  actual  levels  of  heating  in  this 
case  are  about  one  order  of  magnitude  greater  than  in  the  previous  case.  However,  the  position  of  the  maximum 
heating  point  is  the  same  for  both  simulations.  This  should  not  be  a  surprise  since  the  angle  of  attack  is  nearly 
the  same  in  both  cases. 


Comparison  between  DSMC  Results  and  Plume  Heating  Methodology 

During  the  STS- 107  investigation,  a  three-dimensional  plume  heating  methodology  was  developed  based 
on  very  limited  CFD  results  and  represented  a  “highly  engineered”  environment  for  thermal  analysis.  It  was 
desired  to  compare  the  engineering  methodology  to  high-fidelity  DSMC  and  CFD  results  for  STS- 107  types  of 
geometries  and  assess  the  quality  of  the  engineering  predictions  used  for  the  subsequent  thermal  analysis.  Given 
the  complexity  of  the  problem,  the  comparisons  represent  more  of  an  independent  assessment  than  a  validation 
of  the  methodology,  primarily  since  time  did  not  allow  a  second  loop  through  the  process  to  incorporate  lessons 
learned.  Rather,  the  comparisons  focused  on  gross  fluid  dynamic  features  and  qualitative  assessments. 

Two  types  of  comparisons  with  the  DSMC  results  at  350-kft  are  displayed  in  Figure  7.  On  the  left  side  of 
the  figure,  both  sets  of  data  have  been  normalized  by  the  peak  impingement  heating  values  on  the  panel  8/9 
earmuff.  The  DSMC  results  provide  additional  support  for  the  predicted  internal  jet  direction  since  both 
methodologies  predict  peak  heating  values  in  the  same  location.  DSMC  results  also  provide  an  independent 
source  for  secondary  splash  heating  to  the  spar  region  behind  panel  8,  again  in  line  with  the  engineering 
methodology.  Shadowing  of  the  panel  9  spar  region  and  some  enhanced  heating  to  the  panel  9/10  earmuff  are 
also  predicted  by  the  DSMC  results,  in  line  with  engineering  assumptions.  The  right  side  of  the  figure  provides 
a  comparison  of  predicted  heating  magnitudes  with  the  STS- 107  engineering  methodology.  The  engineering 
method  predicts  higher  heating  by  roughly  a  factor  of  two.  However,  the  engineering  methodology  is  based  on 
continuum  assumptions  and  the  calculations  are  made  at  rarefied  conditions  so  the  conservatism  is  not 
surprising. 
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CONCLUSIONS 


DSMC  simulations  of  Columbia  Shuttle  Orbiter  flight  STS- 107  for  the  prime  candidate  accident  scenario 
are  presented.  The  analysis  identifies  instances  where  the  localized  heat  flux  to  the  Orbiter  was  in  the  proximity 
of  the  damaged  leading  edge.  Calculations  of  the  flow  inside  the  leading  edge  also  indicate  the  development  of 
a  recirculation  area  enhancing  the  heat  flux  to  the  walls.  The  DSMC  results  compare  well  with  the  plume 
engineering  methodology  used  in  the  investigation. 


STS  1 07  Engineering  Solution  for  Panel  8  1 0'  hole 


FIGURE  7.  Qualitative  comparison  of  engineering  methodology  with  DSMC  calculations  at  350  kft 
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